#include "Spline.h"

/*void spline::solve()
{
    for (int i = 0; i < order; i++)
    {
        cout << "Points Solved: <" <<  coordinates[i].x << "," << coordinates[i].y << ">"<< endl;
    }
}*/

/*void spline::solve()
{
	double h[order-1], A[order][order], temp[1][1], c[order], b[order-1], d[order-1];
	int  r, co, next=0, check = 1;
                    
            for( r=0; r<(order-1);r++)
            {
                 h[r] = coordinates[r+1].x - coordinates[r].x; 
            }
                 
           for(r=1; r<(order-1); r++)
            {
                 c[r] = 3*(((coordinates[r+1].y-coordinates[r].y)/(coordinates[r+1].x-coordinates[r].x))-((coordinates[r].y-coordinates[r-1].y)/(coordinates[r].x-coordinates[r-1].x)));
            }
            c[0]=c[order-1]=0;
            for(r=0; r<order; r++)                // set all initially to zero
            {
                 A[r][next] = 0;
                 if(r==(order-1) && next != (order-1))      //next column
                 {
                      next = next + 1;
                      r = -1;
                 }        
            }
            A[0][0]=A[order-1][order-1]=1;                 // set 1 Getting ready for Gauss
            for(r=1; r<(order-1); r++)                // assign the values  
            {
                 A[r][r] = 2*(h[r-1]+h[r]);
                 A[r][r-1] = h[r-1];
                 A[r][r+1] = h[r];
            }
            
                                                      
                                                         //Start of Gauss Elimination
      next=1;
      for (r=1; r<3 ; r++)                         // Setting diagonals to 1 and creating upper triangle of matrix
      {
          for(co=3; co>-1; co--)
          {
                      A[r][co] = A[r][co] / A[r][r];
                      if(check==1)
                      {
                                  c[next] = c[next] /A[r][r];
                                  next++;
                                  check = 0;
                      }
          }       
           temp[0][0] = A[r+1][r];
           check = 1;
           
          for(co=1; co<3; co++)
          {
                      A[r+1][co] =  (-1*A[r][co]*temp[0][0]) + A[r+1][co];
                      if (check==1)
                      {
                                   c[co+1] = (-1*c[co]*temp[0][0]) + c[co+1];
                                   check = 0;
                      }
          }
      check = 1;
      }
      
      
      for (r=2; r>0 ; r--)             //  Lower trinagle to complete identity matrix
      {
          temp[0][0] = A[r-1][r];
          for(co=2; co>0; co--)
          {
                      A[r-1][co] = (-1*A[r][co]*temp[0][0]) + A[r-1][co];
                      if (check==1)
                      {
                              c[co-1] = (-1*c[co]*temp[0][0]) + c[co-1]; 
                              check = 0;
                      }          
          }
          check = 1;
      }                                      // after this for loop is finished, the c coeffecients (c[n]) would already-                            
                                             // -be solved
      A[0][0] = A[2][3] = 0;                 // can automatically be set to zero
      
                                            // Note that in form s[r] = coordinate.y + b(x...)+c(x..)+d(....) 
      for(r=0; r<3; r++)                   // b and d coeffecients
      {
                 b[r] = ((coordinates[r+1].y-coordinates[r].y)/h[r]) - (h[r]/3)*(2*c[r]+c[r+1]);
                 d[r] = (c[r+1]-c[r]) / (3*h[r]);
      }
                                              // a coeffecients are automatically assigned to y[n]
      cout<< "The cubic splines are:\n";
      
      for(r=0; r<order-1; r++)
      {
           for(int j = coordinates[r].x; j < coordinates[r+1].x; j++)
            { 
                tempo[r].y = coordinates[r].y + (b[r]*(j-coordinates[r].x)) + pow(c[r]*(j-coordinates[r].x),2) + pow(d[r]*(j-coordinates[r].x),3);
                }
                 cout<<"s"<<r+1<<" = "<<coordinates[r].y<<" + ("<<b[r]<<"(x-"<<coordinates[r].x<<")) + "<<"("<<c[r]<<"(x-"<<coordinates[r].x<<")^2) + ("<<d[r]<<"(x-"<<coordinates[r].x<<")^3)\n"; 
      }
     
      
     
           
}*/

//void spline::graph(double (*func) (double), double LowerBound, double UpperBound, double Iteration)
void spline::graph()
{
    
    for(int r=0; r<order-1; r++)
      {
                cout << "<X,Y>" << r << " is <" << coordinates[r].x << "," << tempo[r].y << ">" << endl;
        }

    
}

void spline::solve()
{
    double A[order-1], B[order-1], C[order-1], D[order-1], temp[1][1], h[order-1];
    double Alpha[order];
    
    int  r;
    //INPUT {(x0, f (x0)), (x1, f (x1)), . . . , (xn, f (xn))}
    //STEP 1 For i = 0, 1, . . . , n - 1 set ai = f (xi ); set hi = xi+1 - xi .
    
    for( r=0; r<(order-1);r++) // getting ready for condition that each of the cubics must join at the knots
            {
                 h[r] = coordinates[r+1].x - coordinates[r].x; 
                 A[r] = coordinates[r].y;
            }
    //STEP 2 For i = 1, 2, . . . , n - 1 set alpha i = 3/hi(ai+1-ai) - 3/hi-1 (ai - ai - 1)
     for( r=1; r<(order-1);r++)
     {
            Alpha[r] = (3*(A[r+1] - A[r])/h[r]) - (3*(A[r]-A[r-1])/h[r-1]);
        }
    //STEP 3 Set Io = 1; set uo = 0; set Zo = 0.
    double I[order], U[order], Z[order];
    I[0] = 1;
    U[0] = 0;
    Z[0] = 0;
    
    //STEP 4 For i = 1, 2, . . . , n - 1 set Ii = 2(Xi+1 - Xi-1) - Hi-1 * Ui-1;
    //set Ui = Hi/Ii; set Zi = (alpha - Hi-1 * Zi-1)/Ii
    for (r = 1; r<(order-1); r++)
    {
        I[r] = 2*(coordinates[r+1].x - coordinates[r-1].x) - h[r-1]*U[r-1];
        U[r] = h[r]/I[r];
        Z[r] = (Alpha[r] - h[r-1] * Z[r-1]) / I[r];
    }
    //STEP 5 Set In = 1; set Cn = 0; set Zn = 0.
    I[order] = 1;
    C[order] = 0;
    Z[order] = 0;
    
    //STEP 6 For j = n - 1, n - 2, . . . , 0 set Cj = Zj - Uj * Cj+1;
    //set Bj = (Aj+1 - Aj)/Hj - Hj*(Cj+1 + 2Cj)/3; set Dj = (Cj+1 - Cj)/3Hj
    for (int j = (order-1); j >= 0; j--) // Solving For Coefficients
    {
        C[j] = Z[j] - U[j]*C[j+1];
        B[j] = (A[j+1] - A[j])/h[j] - h[j]*(C[j+1] + 2*C[j])/3 ; 
        D[j] = (C[j+1] - C[j])/(3*h[j]);
    }
    //STEP 7 For j = 0, 1, . . . , n - 1 OUTPUT aj , bj , cj , dj .
    spline::tempo = new coordinate[order];
    for(r=0; r<order-1; r++)
      {
           for(int j = coordinates[r].x; j < coordinates[r+1].x; j++)
            { 
                tempo[r].y = coordinates[r].y + (B[r]*(j-coordinates[r].x)) + C[r]* pow((j-coordinates[r].x),2) + D[r]*pow((j-coordinates[r].x),3);
                j = coordinates[r].x+coordinates[r+1].x;
                }
                 cout<<"s"<<r+1<<" = "<<coordinates[r].y<<" + ("<<B[r]<<"(x-"<<coordinates[r].x<<")) + "<<"("<<C[r]<<"(x-"<<coordinates[r].x<<")^2) + ("<<D[r]<<"(x-"<<coordinates[r].x<<")^3)\n";
                 
      }
    
}
